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Abstract 

The tilings of lozenges in 2 dimensions and of rhomboedra in 3 dimensions are 
studied when they are constrained by fixed boundary conditions. We establish a 
link between those conditions and free or periodic boundary ones: the entropy is 
written as a functional integral which is treated via a saddle-point method. Then we 
can exhibit the dominant states of the statistical ensemble of tilings and show that 
they can display a strong structural inhomogeneity caused by the boundary. This 
inhomogeneity is responsible for a difference of entropy between the studied fixed 
boundary tilings and free boundary ones. This method uses a representation of tilings 
by membranes embedded in a higher-dimensional hypercubic lattice. It is illustrated 
in the case of 60 degree lozenge tilings. 
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Introduction 



Since the discovery of quasicrystals in 1984 Q, a great deal of work has been accomplished to 
get a precise microscopic structural description of these metallic non-crystalline alloys. It was 
rapidly understood that Penrose-like ||] tilings could provide very good microscopic models 
of quasicrystals: it is highly presumed that favored atomic motifs form tiles. However, the 
best description for real quasicrystals remains an open question: according to the mechanisms 
involved to describe the structure and explain its stability, the studied tilings can be perfect 
Penrose-like arrangements of tiles or random ones. Indeed, despite their random character, 
the latter exhibit global quasiperiodic symmetries || || and are therefore good candidates 
for quasicrystal models. The random tilings use the same tiles as the perfect ones. But 
in the former, local rearrangements of tiles - called localized phason flips - are allowed. 
These degrees of freedom give access to a great number of microscopic configurations. This 
is the random tiling model (RTM) j|, |(| ||. It involves an important contribution of the 
tiling entropy to the total configurational entropy, and therefore to the free energy. This 
phenomenon is supposed to favor the quasicrystal against other competitive phases. 

Among the different techniques developed to estimate this tiling configurational entropy, 
the partition method [7|, ||, ^, [n], 11 1 presents the advantage to set a particularly well defined 



combinatorial problem. However, the boundary conditions of the tilings considered in this 
method are different from the usual ones. As a consequence, the configurational entropy 
per tile of partition tilings is lower that the usual free boundary one: in the simplest case 
of 60 degree lozenge tilings, the respective values can be exactly calculated and are about 



0.261 [0] and 0.323 [12] at the infinite size limit when the 3 fractions of tiles are equal (the 
configurational entropy per tile is simply the logarithm of the number of tilings divided by 
the number of tiles) . 

Elser |?[ explicited a qualitative argument to explain this difference. The goal of the 
present paper is to go further: we will establish the link between the different kinds of 
boundary conditions and we will give a qualitative as well as quantitative explanation for 
the difference of entropy per tile. Hence, we will connect two related models of statistical 
mechanics, which are usually treated via rather unrelated methods. 

A preliminary part of these results were briefly exposed in reference [|10|]. Most of them 
were concisely presented in a shortened version during the Sixth International Conference on 
Quasicrystals (ICQ6), in Tokyo [jnj. 

1 Random Lozenge Tilings and Boundary Conditions 

In this paper, we will consider d-dimensional tilings of rhombic tiles (lozenges in 2D, rhom- 
boedra in 3D) which tile a region of the Euclidean space, without gaps nor overlaps. The 
standard method for generating such structures consists in a selection of sites and tiles in 
a .D-dimensional lattice (D > d) according to certain rules, followed by a projection onto a 
suitable d-dimensional subspace and along a generic direction. We then say that we have 
afl-»d tiling problem. It is the so-called "cut-and-project" method 14, 15]. By con- 
struction, the so-obtained rhombic tiles are the projections of the (f-dimensional facets of the 
D-dimensional hypercubic lattice. There are (j) different species of tiles. In the simple 
3 — > 2 case, this amounts to three kinds of differently oriented 60 degree lozenges. 



2 



Usually, those tilings have periodic or free boundary conditions and it is generally admitted 
that the respective entropies are equal at the thermodynamic limit - for given fractions of 
tiles. As we have just seen it, we will consider another kind of boundary conditions in the 



following, the fixed boundary ones [1C], related to the partition method. Then, the region to 
be tiled will be the generic "shadow" of a D-dimensional rectangular parallelepiped, the sides 
of which take integer lengths in the .D-dimensional hypercubic lattice. This generic shadow is 
called a zonotope |lq] , denoted by Z. When d = 2, the zonotopes coincide with the 2Z)-gons. 
The tiles are supposed to perfectly fit with the boundary dZ of Z. Examples of 3 — > 2 tilings 
are given in figures [l| and |3|. 

Then the entropy per tile is a function of the different fractions of tiles, or in other words a 
function of the side lengths of the zonotope Z. For example, the 3^2 tilings are enumerated 
by MacMahon's formula |jl7|| , which was derived at the beginning of this century: 

, 2 = (fc + f + p-l)!P](fc-l)!P](f-l)!Pl(p-l)!P] 

k ' l ' p {k + P -iy.i 2 ](i+ P -iy.M(k + i-i)\w ' lJ 

Here, this formula has been rewritten in terms of second order generalized factorial func- 
tions ||: 

h 

kl®=k, kl^ = JJi!' m_1 '. (2) 

i=i 

The quantities k, I and p denote the side lengths of the hexagonal boundary. The case when 
all these lengths are equal will be called diagonal in the following, and the corresponding 
entropies will be called diagonal entropies. 

In appendix |A,2 , we rederive this formula via a purely combinatorial method which will 



prove to be useful in the following, the Gessel-Viennot method [18, 

In the following we will be interested in the infinite size limit entropy of such tilings when 
the different fractions of tiles are given. This amounts to making the side lengths of the 
boundary tending to infinity with their relative ratios held fixed: the number of tiles goes to 
infinity while the shape (but not the size) of the boundary is kept fixed. 

2 Continuous Limit and Functional Integral 

2.1 Membrane Representation of Tilings 

This point was already developed in previous publications [l0|, 11] and is closely related to 



the cut-and-project method. Therefore we will only give a brief presentation of this method. 
The main idea is that a random tiling can be lifted as a (f-dimensional non-flat structure 
embedded in a L>-dimensional space. 

This structure is a continuous membrane made of d-dimensional facets of the Z D hyper- 
cubic lattice. When this membrane is projected along the suitable direction, the projections 
of these facets are precisely the tiles the tilings are made of (section |l|); its continuous char- 
acter guarantees the absence of gaps in the so-obtained tiling. Such a membrane is said to 
be directed to emphasize the fact that its projection does not create any overlap. 

For example, figure [l] displays a 3 — > 2 tiling, which can also be seen as a 2-dimensional 
non-flat directed membrane embedded in a cubic lattice. To get a tiling, this membrane 
must be projected along the (1, 1, 1) direction of the cubic lattice. This point of view can be 
generalized to arbitrary-dimension tilings. This correspondence is always one-to-one. 
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Figure 1: 3-dimensional representation of a 3 — > 2 tiling. 

In this D-dimensional point of view, the cf-dimensional space on which the membranes are 
mapped to get the tilings is called the real space £. Its perpendicular space is denoted by £ ± . 
For the sake of convenience, we choose the space coordinates to be the d coordinates on £ and 
the D — d ones on Since the membranes are directed, they can be seen as mono-valuate 
functions <f> from R d to R D ~ d . More precisely, in the case of fixed boundary conditions, these 
functions are defined on the zonotopal region Z of £. 

In the case of free boundary conditions, if the function cj> has a large scale global gradient 
V4> = E, the random tiling model states that the fractions of tiles are controlled by this 
gradient and, therefore, the entropy per tile can be written as a function of E Q. This 
gradient is usually called the phason gradient. 



Boundary Conditions in the Membrane Representation 

We now have to define what the fixed boundary conditions of section [l] become in the 
language of directed membranes. As illustrated on figure |l|, we also get a boundary condition 
in the -D-dimensional space: the membrane is inscribed inside a (non-flat) polygon (or poly- 
hedron), the projection of which on the d-dimensional space gives the tiling boundary, dZ 
TO , Oj. For instance, the boundaries of 3 — > 2 tilings are (non-flat) 3-dimensional hexagons. 
We will call these boundaries the membrane frames. Such a frame can be precisely charac- 
terized as the inverse image, via the projection, of the boundary of the zonotopal region Z 
which is being tiled [11]. It is therefore a subset of the boundary of the hypercube from which 



Z originates. This frame results in conditions on the functions (f>, on the boundary dZ of Z. 

In the case of free or periodic boundary conditions, the functions (j) have free or suitable 
periodic conditions on the boundary of the domains T> on which they are defined (here, we 
use the notation D instead of Z to denote domains which might be non-zonotopal) . 



2.2 Continuous limit 

Once a tiling has been "lifted" in the higher-dimensional space, the so-obtained directed 
membrane has a corrugated aspect, due to its discrete character. Here, we wish to get rid of 
this discrete character, at the infinite size limit, in order to study more regular ("smoother") 
objects, to which analytic tools can be applied. Moreover, these objects will turn out to 
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characterize the macroscopic states of the statistical ensemble of tilings (or membranes). Let 
us explain how this continuous limit is taken. 

So far, we have considered tiles of side length 1. To define the continuous limit, we will 
get rid of this discrete character, thanks to a suitable rescaling. For reasons that will become 
clear in the following, this side length will go to as the number of tiles tends to infinity. 
The functions which represent the tilings are defined on a domain T> of the real space. If TV" 
is the number of tiles in T>, we need to rescale the tilings by a factor \/N l / d . Thus in any 
infinitesimal domain d d y in T>, the number of tiles goes to infinity when N does. Moreover, 
we do the same rescaling in the perpendicular space £ . 

Once we have done this rescaling, the tilings are represented by functions </> which have 
quite an irregular aspect at small scales. As it is usually done, for instance in polymer or 
polymerized membrane theories, we will treat small scale fluctuations and large scale ones 
in a different way. Large scale fluctuations are represented by regular functions whereas 
microscopic ones around the latter functions are integrated in an entropic term. The latter 
term will have an exponential form, and will therefore be treated via classical methods on 
functional integrals. 

Below, we will adopt the following terminology: a membrane (or function) which has 
microscopic fluctuations, that is to say which is the exact representation of a large tiling, will 
be called faceted. A membrane, the fluctuations of which have been integrated in an entropic 
term, will be called smooth. Finally, we will go on calling a tile any d-dimensional facet of a 
faceted membrane. 

Given a smooth membrane 0, we can adopt the first nave definition of the entropy of this 



membrane [10]: 



, , log(Number of iV-tile faceted membranes close to 0) 

sM = & n • (3) 

The important point here is to understand that, thanks to the previous rescaling, (j> is 
kept fixed while the number of tiles goes to infinity. This point allows, among many other 
things, to work with the same set of functions whatever the system size. 

In this definition, we have not characterized what "close" meant. Let us first state that 
its precise definition is unessential. To understand this point, let us focus on an analogy with 
a far more simple statistical mechanics problem, where key ideas are easier to catch: let us 
consider a closed box containing a perfect gas. This box is divided into two parts of same 
volume, A and B, separated by a virtual frontier. Then one looks for the entropy &{xq) of 
this system, at the thermodynamic limit, when a fraction xq of the molecules lie in A (and 
therefore a fraction 1 — xq in B). Of course, this quantity xq fluctuates and can only be 
defined up to a small quantity Ax. Then we write 

log (Number of configurations such that x = xq db Ax) 
a(x ) = hm , (4) 

where N is the number of molecules. The important point here is that it can then be proven 
that this definition of cr(xo) does not depend, at the thermodynamic limit, on the precise 



choice of Ax, provided it is a finite quantity (see [ 20 ] , p. 30, for example, for a discussion on 
this point). 

In particular, if one looks for the more likely value of xo, that is for the maximum of cr(xo), 
one finds xo = 0.5, still independently of the choice of Ax. In other words, the dominant 
states are such that x is close to xq, but is not exactly equal to xq. 
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This point of view also applies to directed membranes: we will see in the next paragraph 
that this is all the more justified since at the infinite size limit, only a few constraints cause 
the faceted membranes to be "stuck" to the smooth one, 4>. 

The next step consists in considering a point yo and an infinitesimal domain d d y in T> 
containing yo- This domain is large as compared to the tile size, since this size tends to 0. 
Moreover, since d d y is infinitesimal, the gradient V</> can be considered as constant on this 
domain. Therefore d d y contains a piece of tiling with an "infinite" number of tiles and a fixed 
phason gradient E = V</>(yo)- This phason gradient is the only constraint on this piece of 
tiling. In particular, its boundary conditions are free. Hence if <r{E) denotes the entropy per 
tile of a large free boundary tiling of global phason gradient E, then the number of faceted 
membranes close to <f> and defined on d d y is equal to 



M(y ) = exp iV^MVcKyo)) 



(5) 



where N{d d y) is the (large) number of tiles of the previous membranes. This number of tiles 
depends on the domain size (and on the total number of tiles, N): 

N(d d y) = N.n{Vcp)d d y, 

where n(V0) is a function, the integral on T> of which is equal to 1. Hence 



M(y ) = exp N.n(V<f>(y ))(T(V(t>(yQ))d d y 



Hence the total number of membranes close to on the whole domain D is given by 

N<t> = n^(y) = II ex P k^v<A(y)MV0(y))d d y 



(6) 



(7) 



This product runs on all the infinitesimal domains d d y. Rigorously, since the membranes 
are to coincide on the boundaries between the different domains, this product should be 
divided by a boundary correction term. But when N —* oo, the latter infinitesimal domains 
contain an infinite number of tiles and these boundary terms disappear]]. The total number 
of membranes is therefore equal to the product of the individual numbers of membranes in 
each domain d d y when N is large. 

In order to simplify the forthcoming calculations, we choose the infinitesimal domains so 
that they form an hypercubic lattice which divides the large domain T> in small cubes of side 
a and of volume a d = d d y. These domains are indexed by Greek letters and their vertices by 
Latin letters. 

We have already seen that if a is small enough, then <j> can be considered as afhne on any 
domain a. Therefore this function only depends on its values (f) a j on the vertices % of the 
domain: 



Ms 



J] exp [iV.n(£(<^)a(£(<^)d d y 

a 

Y, [N.n(E((t>a,i))<r(E(<f> aA ))d d y 



exp 



1 In other words, the entropy of two infinite size sub-systems in which V<j> is fixed is additive. 
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which is actually a function M{(pi) of the values fa on the lattice vertices. 

To sum up, we have fixed the global shape of the faceted membranes: they are tied to 
the lattice vertices, that is to say we have imposed their mean gradient on the domains a to 
be equal to the gradient of 4>. It is actually the only constraint we have imposed to these 
membranes which are counted by Ms- Now, the key point is that this constraint is sufficient 
to be sure that, at the infinite size limit, the faceted membranes counted by Ma, tend towards 
cf>. Indeed, Henley Q showed that for any dimension d, if the gradient E is fixed, if Ah(L) 
denotes the fluctuations in the perpendicular space of directed membranes of linear size L, 
then 

Ah{L) 



L 



when L — > oo. 



L 1 / 2 , if d = 2, then Ah(L) 



(8) 

log L and if d = 3, 



(More precisely Q, if d = 1, then Ah(L) 
then Ah(L) is uniformly bounded). 

After the l/N l / d ~ 1/L rescaling, these fluctuations tend to when N — > oo. Therefore 
all the membranes counted by M^ tend uniformly towards (j). Hence, they are close to </>, 
whatever the precise definition of this term. Finally, 



lim 

N— >oo 



log(AQ 



N 



Now, the total number of faceted membranes, M, is given by the integral^]: 



M 



IIWi)= / Hd<p i exp^ N [ n ( E ^c e , l )ME( ( / )a:i ))d c 



0) 



(10) 



So far, we have discretized the domain T> to be sure that the membranes have an infinite 
number of tiles in any infinitesimal domain a, at the infinite size limit. Then A/^ or J\f count 
faceted membranes close to smooth membranes which are afEne on every such domain. To 
get rid of this restriction, we will now take the limit a — > 0. Formally, we write 



Vd> = lim TT 



and we turn the sum into an integral. Thus 



and 



Therefore 



P^exp 



N / n{V(/))a(V<p)d d y 
v 



Md> = exp 



N I n{V(/))a(V<p)d a y 
v 



n(y<j>)a(V<p)d d y. 

v 



(11) 
(12) 
(13) 



As foreseen, this expression is independent of the precise characterization of the above "close- 
ness" in the definition of s\d>]. 



Rigorously, in the case of free or periodic boundary conditions, the problem is invariant under translations 
of 4>i an d this integral M is divergent. The membrane must be fixed to a point to avoid this divergence. For 
example, we fix <j>(y = 0) = 0. 
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This way of writing the entropy associated with a smooth membrane <p is quite similar 
to Henley's (section 6.1). He gets this result thanks to a suitable coarse-graining of 
the membranes. The coarse-graining of a faceted function is essentially its local mean in a 
neighborhood of diameter ao of every point of £>(Q). Nonetheless, in his point of view, ao 
is large but finite, whereas the smooth functions that we consider here integrate the local 
fluctuations of an infinite number of tiles when N becomes infinite. Moreover, there is a 
technical difference in the two expressions: Henley's entropy is an entropy per unit area, 
whereas ours is an entropy per tile. The ratio between these two entropies is actually n(V(f>), 
the tile density (per unit area). 

Indeed, this last quantity directly depends on the phason gradient E = V</>, since the 
latter controls the different fractions of tiles and the different tiles do not necessarily have the 
same area. The entropy per unit area, n(V</>)<r(V0), will be denoted by t(V0). 

However, in the codimension-one case, all the tiles have the same area or volume, and 



equation 13 can be simplified. The tile density, n, does not depend on the phason gradient 
any longer and can be factorized before the integral. The exact value of n depends on the 
choice of the rescaling: so far, we have only specified its order of magnitude (l/N 1 ^) but 
not its exact value. To calculate n, we choose <j) to be zero everywhere. Then = and 
cr(V^) = o"o- The entropy per tile s[4>] is also equal to ctq, since this membrane has a vanishing 
gradient. Hence o"o = n j^aod^y and n = 1/V(T>), the inverse volume of T>. 
In the codimension-one case, 



/ a(V0)rf 
Jv 



"y 



*1 = viv) ■ (14) 



which is the expression we had given in reference pi 



To sum up, we have coded the macroscopic states of this statistical ensemble by an 
internal parameter 4>, and we have calculated the entropy associated with these states. Let us 
emphasize that the functional s[(f>] is expressed in terms of the free boundary tiling entropy 
a, whatever the conditions on the boundary dT> of the domain T>. An example of smooth 
function <f> is displayed in figure ||. 

2.3 Dominant States in the Statistical Ensemble 

The faceted membranes defined on T>, and that have the good boundary conditions, are 
counted by 

M = [ V<\> exp[Ns[4>}}. (15) 

Let us precise that this functional integral is taken upon the set F of functions which are 
smooth representations of faceted membranes. 

Moreover, the entropy per tile associated with this set of all the membranes (or tilings) is 
given by 

S=lim*f^. (16) 

3 More precisely, the coarse-graining is a convolution product of (j> and another function of spatial exten- 
sion ao. 
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Figure 2: A smooth function (f> : R 2 — > and its frame, in the 3 — > 2 case. 

Now, we suppose that there exists an unique function </> max that maximizes the entropy 
functional s [</>]. This point will be discussed at the end of this section. Moreover, we suppose 
that the functional is non-singular near this maximum, that is to say it has a quadratic 
behavior: 

S[4>] = S[(f> max ] - d d U d d V k u , v (<p(u) - max (w), <p{v) - max (u)) + . . . , (17) 

Jv Jv 

where k is a positive quadratic form, which a priori depends on the point (u, v). 

Hence, near <^ max , AA is a Gaussian functional integral, and thanks to a generalized saddle- 
point argument, we get 

8 - lim 

N^oo N 
= S [</>max] • 

This is a classical result in statistical physics: at the infinite size limit, the total entropy is 
equal to the entropy of a dominant macroscopic state. In other words, the statistical ensemble 
of faceted membranes is dominated by membranes close to (p max . In the space of membranes, 
the distribution is "peaked" around </> max at the infinite size limit, and looks more and more 
like a Dirac distribution. 

To close this section, we must discuss the assumption of unicity of max - We will use a 
general convexity argument: if a function / is strictly concave on a convex set C and if / has 
a maximum on C, then this maximum is unique. 

Now, the set F of functions is convex: whatever their boundary conditions on &D, let 4>i 
and 4>2 be any two smooth functions in F and let A be any real number between and 1, then 
4>\ = X(f>x + (1 — \)4>2 is also an element of F. In particular, if Vcfti and V<fo satisfy the good 
conditions to insure <j>\ and <p2 to be in F, by linearity of the gradient, V(A</>i + (1 — X}^) 
satisfies the same conditions. And whatever the boundary conditions, (f>\ also satisfied them. 

Let us now check whether s[<j>] is concave: it is an integral, and therefore a positive linear 
combination of functions of <j> (the entropies per unit area) . If we prove that all these functions 
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are concave, then it follows that s[(ft] is concave. Now, the entropy at the point y in T> is the 
composition of t(E) and the function (ft V^>(y), which is in turn a linear function. Hence 
we only need to prove that t{E) is a concave function of E. 

This property is a little stronger than the general random tiling model hypothesis, which 
states that the free entropy has a unique maximum as a function of the gradient E, and is 
quadratic near this maximum Q: 

T f ree {E) ~ r /ree - ^Rf ree E 2 , (18) 

where K^ ree is the so-called tensor of phason elastic constants. Even if our stronger hypothesis 
is a priori valid for a more restricted set of models, let us note that it is satisfied in all exactly 



solvable tiling models [21, 22, p3[ |. The concavity of s[<f>] and the unicity of (ft max are therefore 
reasonable hypotheses. 

Finally, let us remark that since this maximum is unique, it will respect all the underlying 
symmetries of the problem. We will see an illustration of this fact in the following. 



3 Relationship Between Different Boundary Conditions 

In principle, whenever the free boundary entropy r is known, the functional s[(ft] is precisely 
defined, and one can therefore get the fixed boundary entropy. Theoretically, we are able 
de deduce the maximum entropy of fixed boundary tilings, TQ txed , as well as the phason 
elastic constants, K^ lxed , from their counterpart in the free boundary case, Tq 7 * 66 et K^ ree , 
and to invert these relations. This was done in the 3 — > 2 case, in reference J0|(Q), in the 
so-called quadratic approximation, which consists in estimating the free boundary entropy by 
its quadratic development (equation |l^) near its maximum. Here, we will go further and give 
a complete treatment of this 3^2 case. We will also present some related work in the case 
of different fixed boundary conditions. 

To go beyond the quadratic approximation, we will characterize the maximum of the 
entropy functional s[cft] by means of a functional derivation. If s[(ft] = fj ) n(V(ft)a(V(ft)d d y, 
then 



V(V) 
Hence, 



Ss = s[(p + 5(j)} - s[4>] 

l — j^dT(V^(y)).V(5^d d y (19) 

6<ft div{dT)d d y. 



V(V) 
-1 



v 



-div(dr(V<Ky))). (20) 



s<Kv) 

Therefore <^ max is the function (ft which satisfies this equation and which has the good boundary 
conditions. 



4 In this reference, the value 0.253 of the diagonal entropy in the quadratic approximation was erroneous, 
because of badly controlled boundary effects. The actual value is 0.251. 
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3.1 Hexagonal Tilings 

Expression is general. In the 3 — > 2 case, it reads 

d 2 r d 2 (f> d 2 r d 2 <j) d 2 r d 2 cf> , , 

_ -l 9 — I — = f) (21) 

8E 2 dx 2 8E 1 dE 2 dxdy 8E 2 dy 2 1 ' 

where E = {E 1 ,E 2 ). 

d 2 r d 2 r d 2 r 
The coefficients — — ^ , — ■ — ■ and — — k are known. Indeed, in the 3 — > 2 case, there exists 
dE( 8E 1 dE 2 dE^ 

an analogy between the tilings and the ground states of an antiferromagnetic Ising model on 



a triangular lattice [21|. The entropy can then be derived from the previous solution of this 
Ising model [24]. In the latter reference, the entropy is written in terms of chemical potentials. 
Some algebraic manipulations are therefore necessary to write it in terms of E, and then get 
the above coefficients p| : 

8 2 t vr 1 fl + w 2 \ 

77-^7 = 77 T-COS0 , (22) 

8E 2 9sm6\l-w 2 J 1 V ; 

d 2 r _k w_ 2-w 2 

dE 1 dE 2 3^3 sin9 1 - w 2 ' { ' 

d 2 r 2ir w 1 

dE\ 3 1 — w 2 sin <j) ' 

where 6 = ^-^ L (E 1 + y/2), 4> = J\kE 2 , and w = tan(6»/2)cotan(0/2). 
o V <j 



(24) 



The partial differential equation 21 can be solved by means of numerical calculations. The 
idea is to discretize the domain Z, which is a hexagon in this case, and to use an iterative 
process: at each step, a function (f>k is computed. The above coefficients are calculated in 
terms of Then (pk+i is the solution of a linear system which is the discrete version of 



equation 21. Then the sought function is the fixed point of this iterative process, which is 
reached after about 10 iterations. This method was used in reference jHJ and gave satisfactory 
numerical results. 

However, more recently, we were aware of related works by mathematicians which are 
interested in similar problems. They are indeed able to exactly compute the function max 
by the mean of purely combinatoric methods 2S, p6| : the idea is to calculate the number 



of fixed boundary tilings with a precise distribution of vertical tiles upon a given horizontal 
line. The value of m ax on this horizontal line is then given by the distribution of vertical 
tiles which maximizes the number of such tilings at the infinite size limit. 

This exact solution points up a very singular phenomenon. The above authors called 
it the arctic circle phenomenon: at the infinite size limit, in the tiling representation in 2 
dimensions, there is a central region of £ which is circular in the diagonal case and elliptic in 
general, inside which the tiling is random - in this sense that it contains the 3 kinds of tiles. 
Outside this region, the tiling is "frozen" : as illustrated in figure ||, there are 6 regions where 
there are only one kind of tiles, and where the entropy is equal to 0. 

In these regions, the function (/> max is rigorously linear and its gradient is constant. Inside 
the central region, in the diagonal case, the phason gradient of ^ max is equal to 



E x = -V2 + —= [cotan- 1 f(x, y) + cotan- 1 /(-x, y)] , (25) 

7TV 2 
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Figure 3: A randomly generated 3^2 tiling. At the infinite size limit, there is a circular 
central region (the "arctic circle" |25[| ), where the tiling contains 3 kinds of tiles, and 6 
"triangular" regions, where the tiling contains only one kind of tiles and is said to be "frozen" . 




E 2 = -\ - [cotan 1 f(x, y) - cotan 1 f(-x,y)] , (26) 



where f(x,y) = — = ^Z_ J^ x ^ ^Zj ^ ~*~ if the side of the hexagon has been rescaled to 

Jy ,y) 2^3 Vi- 4 /3(^ 2 + y 2 ) 

1. In this expression and the previous ones, the origin of the coordinates is at the center of 
the hexagon, the axis (Ox) is pointed towards a vertex of the hexagon. 

These expressions characterize ^> m ax- In fact, the resulting function is very close to the 



function showed in figure 0, which was actually computed in the quadratic approximation 17 



[10]. Now, it is possible to compute the entropy per tile associated with max : a simple 
numerical calculation gives S[0 max ] = 0.2616(3). This value is in exact agreement with 
the exactly known diagonal entropy of fixed boundary tilings. These results are therefore a 
validation of our continuous approach (coarse graining). As announced, we have established 
an exact link between free and fixed boundary tilings. 

This function </> max deserves a quick qualitative description; since ^> max is very close to the 
function displayed in figure ^, we will use this figure to illustrate our arguments: first, when 
the boundary is an hexagon, it has a 3-fold symmetry and, as foreseen, the solution </>niax 
respects this symmetry. Second, because of the strong influence of the boundary, there is a 
gradient of entropy between the boundary and the bulk. Indeed, near the center of the tiling, 
the gradient of <^ max is very close to the free boundary one, whereas far from the center, this 
gradient becomes more and more influenced by the boundary and becomes singular out of 
the arctic circle: there, the entropy is zero. Then the fixed boundary tilings provide a very 
interesting model having an inhomogeneous entropy distribution. 

To close this discussion, let us precise that this infinite size limit cannot be called a 
"thermodynamic limit" because of this lack of homogeneity: in statistical mechanics, a system 
is said to be at the thermodynamic limit if its properties do not depend on how it tends to 
infinity. In particular, they must be homogeneous in the system and must not depend on the 
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container shape (here, the boundary) j20|. Here, the situation is far from this one, since even 
the stoichiometry depends on the boundary shape. 



3.2 Other Kinds of Tilings 

So far, we have only considered zonotopal fixed boundaries. The reason for this choice is that 
the first motivation of this work was to establish the link between free boundary tilings and 
tilings built by means of the so-called partition method, the aim of which was to develop a new 
approach to these tiling questions (7|, |8|, |9|, 11 ] . The latter tilings precisely have zonotopal 



boundary conditions, by construction. However, there is no reason why the previous method 
could not be applied to other kinds of boundaries. In this section, we will analyze tilings, the 
boundary of which is fixed, but which nonetheless have a free boundary entropy. 

Indeed, if the boundary is fixed but imposes a uniform phason gradient, that is to say the 
membrane frame lies in a d-dimensional plane of gradient E = Eq, then the affine membrane 
of gradient Eq has on the one hand the good boundary conditions, and on the other hand 
satisfies the partial differential equation ^l] (since its second order derivatives are equal to 
zero). This membrane is therefore the function (/) max and its entropy s[</>max 

] is equal to the 

free boundary one, t(Eq). 

More precisely, let us analyze a 3 — > 2 class of tilings, the fixed boundary of which is Hat 
in the membrane representation in 3 dimensions, that is to say does not impose any global 
phason strain to the tilings. This kind of boundary is illustrated in figure ^. 



Figure 4: Left: a zigzag boundary, without phason strain. Right: tiling detail; the boundary 
(thick line) is globally - but not locally - strait. 

At the infinite size limit and in the membrane point of view, after rescaling, the boundary 
becomes a flat hexagon, that is the function (f> is constrained to zero on this boundary. Then 
the function eft = maximizes the entropy and, at the infinite size limit, this entropy is equal 

to 4 ree . 

We have tested this theoretical prediction by a direct calculation of this entropy. The 



method is developed in appendix |A.3j . It uses again the Gessel-Viennot method [18, IS]. The 
entropy of very large tilings can then be numerically reached and fitted to get its infinite size 
limit. We find an entropy per tile of 0.32306(4), which is precisely the free boundary entropy. 

To close this section, let us draw attention to numerical simulations by Joseph and 



Baake [27]: they analyzed the configurational entropy of random 4 — ► 2 tilings, the boundary 
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of which is fixed and flat (in 4 dimensions), as in our previous example (the global phason 
strain is zero). As foreseen, the entropy that they eventually found was equal to the free 
boundary one, which was itself numerically estimated. 



Conclusion 

Thanks to a continuous approach in the membrane point of view at the infinite size limit, 
we have been able to describe the dominant states of statistical ensembles of tilings. In 
particular, this method has enabled us to establish a quantitative link between two exactly 
solvable models of statistical mechanics, the free and fixed boundary tilings of 60 degree 
lozenges. In the latter, a very remarkable event occurs, the arctic circle phenomenon: there 
exist "frozen" regions of the tilings in which there is only one kind of tiles and were the entropy 
is therefore zero. This lack of homogeneity is responsible for the difference of entropy between 
the two problems, even if they were at first sight closely related. Moreover, this infinite size 
limit cannot be called a thermodynamic limit because of this lack of homogeneity. 

In the case of larger dimension or codimension tilings, a similar treatment would require 
the knowledge of the free boundary entropies. Unfortunately, despite a great deal of work, 
these entropies are not known yet. However, the maximum (diagonal) entropies and the 
phason elastic constants are numerically known in many cases. It could therefore be possible, 
in these cases, to compute the fixed boundary entropies and phason elastic constants, in 
the quadratic approximation. But this approximate method would not be conclusive on the 
existence of an arctic circle phenomenon in such problems, which is nonetheless a captivating 
open question. 

Finally, it is worth emphasizing that this method could be useful in describing how any 
constraint imposed at the boundary relaxes into the bulk. In this paper, we have studied 
two kinds of boundaries. The first one, the straight boundary related to partition problems, 
imposes to the tiling the strongest constraint among all boundary conditions: the tilings must 
relax continuously from a completely crystalline structure to a random one. Physically, such 
tilings (in 3D) could model the result of a growth of quasicrystalline materials on crystalline 
phases. More generally, more complex physical situations, such as extended topological defects 
(such as dislocations), or other kinds of interfaces, could a priori be translated in suitable 
boundary conditions. The numerical method we have developed could then be applied to any 
such boundary conditions and could be useful in describing how the material relaxes in the 
presence of such constraints. 



Note: we have recently been aware of related works in the Aztec diamond tiling problem [28]. 



A The Gessel-Viennot Method 

In this appendix, we present a combinatorial method for the counting of configurations of 
avoiding paths on planar graphs, the Gessel-Viennot method, which can be very useful in the 
enumeration of fixed boundary tilings. We illustrate this method in two 3 — > 2 cases discussed 
in the present paper. 
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A.l The Method 



We will not extensively present the Gessel-Viennot method [18, 19]. We will instead give a 
brief description and try to explicit the underlying ideas. The method is rather general: it 
can be applied to any oriented graph without cycles (acyclic oriented graph), in which are 
selected two families of vertices, Ui and V{, i = 1, . . . , n. This graph is supposed to satisfy the 
property of compatibility: if two paths on this graph are going respectively from to Vj ± 
and from Ui 2 to Vj 2 and if these paths do not cross, then i\ < 12 and j\ < ji- Note that this 
property is very specific to two-dimensional graphs. 

We are interested in the number of configurations of n avoiding (or non-intersecting) 
paths, the i-th path going from Ui to V{. If we denote by Ay the number of paths going from 
Ui to Vj, then the method states that the number of configurations is equal to the following 
determinant: 



det(Ajj)i<j ) j< n . 



(27) 



The idea of the proof is the following: in this determinant, all path configurations, whether 
intersecting or not, the i-th path going from U{ to t> CT (j), for any permutation cr, are counted, 
with a + or — sign. All intersecting configurations cancel two by two and only the non- 
intersecting configurations remain. They are exactly the sought configurations thanks to the 
property of compatibility. The reader interested in more details can refer to the review paper 



by Stembridge ]19|. 



A. 2 Hexagonal Boundaries 

In this section, the previous method will be used to rederive MacMahon's formula (see sec- 
tion [l]). Consider a 3 — > 2 tiling of an hexagonal region of side lengths k,l and p (figure 
left). In such a tiling, one can follow sequences of tiles which have a horizontal edge. These 
lines, which are usually called worms, cross the tiling from bottom to top. The tiling can now 
be slightly deformed so that a kind of tiles become squares (figure ||) and the p worms can be 
seen as up-going paths on a square lattice (figure [5L right). 




Figure 5: A hexagonal boundary 3 — > 2 tiling (left): the dashed worms can be translated into 
a set of p avoiding oriented paths on a square lattice (right). The i-th path starts from the 
(fixed) vertex Ui and goes to the (fixed) vertex V{. There are p paths. The side lengths of the 
hexagon are p, k and I. 



Therefore the previous theory on avoiding paths on an acyclic oriented graph can be 
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applied. The number \j of paths joining the vertices Ui and Vj is a binomial coefficient: 

(k + l)\ 
(k + j-i)\(l + i-j)V 

Then, we have to compute the following determinant: 

(k + iy. 



D p (k,l) = det(Ajj) = det 



l(k + j-i)\ + 



l<i,j<p 



(28) 



(29) 



D p (k,l) = [{k + l)\] p det 

= l(k + iW 

x det 



1 



(k + j-i)\ + 

1 1 



l<i ,j<p 



(k + p-l)l...k\ (l + p—l)\...U 
{k + p-i)\ (i + i-l)!" 



l<i,j<p 



(30) 



(k - iy.w a - iy.w 

The first factor equals \(k + i)!P — —r-r — where we have used again 

(k + p — ly.w (I + p — 1)H J 
the second order factorial function. The second factor is denoted by A p . We will use the 



notation: (i) 



{k+p-i)\ (l + i-iy. 



(p) 

. Since j < p, P. is a polynomial of degree 



(P ~ j) + (j ~ 1) = P ~ 1- 

We now use the following result concerning polynomials: if Qj, j = 1, . . . ,p are polyno- 

mials of degrees smaller than p — 1, if Qj = ^ a\ X 1 , and if real numbers, 



then 



i=0 



det[Qj(x i )] 1 <i t j<p = det(a| i) ) x VdM(xi, . . . ,x p ), 



(31) 



where VdM(xi, . . . ,x p ) is the Van der Monde determinant of these real numbers. We recall 
that 



VdM(xi, . .. ,x p ) 



1 2 p— 1 

1 Xl X\ ... x x 



1 X2 X 



1/ v* 'V' 2 



(32) 



The proof of equation 31 is straightforward: the left-hand side matrix is the product of the 
coefficient matrix and of the Van der Monde matrix. Note that VdM(l, 2, . . . ,p) = (p — I)!' 2 ! . 

Now, the calculation of A p is made by induction on p: if d p denotes the determinant 
of the coefficients of the polynomials Pj P \ j = 1,... ,p, then thanks to the above results, 
A p+ i = pl^d p+ i. After a tedious calculation, one gets that 



x p+i 



So by induction on p, 



P&Pik+p + l) = (k + l+p)...(k + l + l). 



+ l+p-l)]W 



[(k + iy.] p (fc + j-ijira 



(33) 



(34) 
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(k + i + p-iy.w (k-iy.w (i-iy.w (p-i)g 
pl ' j " (fc+ P -i)!P] (i+p-iym ( fc + /-i)![2] ' ^ 



Finally, we get 



which is precisely Mac Mahon's enumerative formula |17|| , rewritten in terms of generalized 
factorials. 

A. 3 Flat Fixed Boundaries 

This method can also be applied to the tilings studied in section |3T2| (figure |||). However, 
in this case, we will not be able to get an explicit enumerative formula and some numerical 
calculations will be necessary to have access to the infinite size entropy. 

With these boundary conditions, the worm and path representations of figure [5] must be 
modified, as illustrated in figure |(| 




Figure 6: The worm representation of a flat-boundary tiling (left) and its counterpart in the 
avoiding path representation (right), in the case when the number of lines n is equal to 5. 



The two main differences are the following: first, the vertices U{ and V{ are not as simply 
distributed as in figure |5[ Second, some of the n paths are constrained by the presence of 
two vertical bounds (thick dashed lines). Therefore, the number \j of paths starting from 
Ui and going to Uj might be different from a simple binomial coefficient. This number can 
be calculated thanks to the usual "mirror" or "image method" (see [29j, for instance). In the 
diagonal case, with the indexation of figure and when n is odd, 



A; 



Pij 



Pij 

Pij+3(i+j)-2 



Pij 

Kj +6n-3(i+j)+4 
2 



where 



Pij = 2n 



n + 1 



3 



n + 1 



(36) 



(37) 



is the length of any path going from U{ to Vj . 

The method is now strictly similar to the previous one. However, the complete calculation 
of the determinant det(Ajj) seems out of reach. That is why we have chosen to compute it 
numerically for large systems. The so-obtained values are displayed in table |]. In this table, 
n still denotes the number of worms. 
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Table 1: Entropy per tile of n-worm tilings with flat boundaries. These boundaries do not 
impose any global phason strain. 



Number of worms 



Entropy per tile 



n = 21 
n = 31 
n = 61 
n = 81 
n = 101 
n = 151 



S = 0.311881 
S = 0.315379 
S = 0.319098 
S = 0.320065 
S = 0.320653 
S = 0.321446 



The last four data can be fitted with the following law: 




(38) 



Then we get So = 0.32306(4), which is the infinite size entropy (and A ~ —0.246, B ~ 0.245). 
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